Accepted for publication in ApJ Supplement 

Preprint typeset using L^T^X style cmulatcapj v. 6/22/04 



O 
O 

(N 



in 

(N 



MODELS FOR NONTHERMAL PHOTON SPECTRA 
John C. Houck and Glenn E. Allen 

MIT, Kavli Institute for Astrophysics and Space Research, 77 Massachusetts Avenue, Cambridge, MA 02139-4307 

Accepted for publication in ApJ Supplement 

ABSTRACT 

We describe models of nonthermal photon emission from a homogeneous distribution of relativistic 
electrons and protons. Contributions from the synchrotron, inverse Compton, nonthermal brems- 
strahlung and neutral-pion decay processes are computed separately using a common parameteriza- 
tion of the underlying distribution of nonthermal particles. The models are intended for use in fitting 
spectra from multi-wavelength observations and are designed to be accurate and efficient. Although 
our applications have focused on Galactic supernova remnants, the software is modular, making it 
straightforward to customize for different applications. In particular, the shapes of the particle distri- 
bution functions and the shape of the seed photon spectrum used by the inverse Compton model are 
defined in separate modules and may be customized for specific applications. We assess the accuracy 
of these models by using a recurrence relation and by comparing them with analytic results and with 
previous numerical work by other authors. 

Subject headings: radiation mechanisms: non-thermal — supernova remnants 
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1. INTRODUCTION 

Cosmic-rays with energies up to about 1000 TeV are 
thought to be accelerated in supernova remnant shocks 
(Gaisser 1994). In the simplest, idealized model of the 
acceleration process, test particles interact with a shock 
discontinuity to produce a power-law momentum dis- 
tribution of cosmic-rays. Photon spectra produced by 
power-law momentum spectra, as in the XSPEC models 
SRCUT and SRESC (Reynolds 1998; Reynolds & Keohane 
1999), have been widely used to fit the X-ray synchrotron 
spectra of several SNRs. 

However, in the diffusive shock acceleration picture 
(Blandford & Eichler 1987; Bell 1987; Ellison & Reynolds 
1991; Berezhko & Ellison 1999), nonlinear processes are 
expected to cause some deviation, or curvature, away 
from a power-law momentum distribution. Recent obser- 
vations of Cas A (Jones et al. 2003) and SN1006 (Allen, 
Houck & Sturner 2005) have shown evidence for a curved 
synchrotron spectrum. Such spectra can be studied us- 
ing simulations of diffusive shock acceleration in super- 
nova remnants (see e.g. Ellison & Reynolds 1991; Baring 
et al. 1999; Bykov et al. 2000). These simulations predict 
the momentum distribution of nonthermal particles and 
the resulting photon emission spectrum. Unfortunately, 
such detailed calculations are still too time-consuming for 
widespread use in iterative fitting of models to observed 
spectra. 

Therefore, we have developed models of the syn- 
chrotron (§3), inverse Compton (§4), nonthermal brems- 
strahlung (§5) and neutral-pion decay (§6) spectra pro- 
duced by homogeneous emitting regions having non- 
thermal particle momentum distributions with arbitrary 
shape. Our primary goals in developing these models 
were to make them as general and accurate as possible, 
and to make them computationally efficient enough to be 
practical for use in iterative fitting of observational data. 
An additional goal was to adhere to a modular design 
to simplify customizing the models for specific applica- 
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tions. By providing alternate implementations of the ap- 
propriate modules, one can customize the shapes of the 
particle distributions and the shape of the photon spec- 
trum used to compute the inverse Compton model. For 
example, Allen, Houck & Sturner (2005) used the syn- 
chrotron spectrum model with the curved particle mo- 
mentum spectrum described in §2 to detect curvature in 
the cosmic-ray electron spectrum of SN1006. 

We assess the accuracy of the computed photon spec- 
tra by applying a recurrence relation and by comparing 
with analytic results and with published numerical mod- 
els (§7). When using these models to fit (§8) simulta- 
neously radio, x-ray and gamma-ray observations of su- 
pernova remnants, we have occasionally found it useful 
to reduce the number of degrees of freedom by introduc- 
ing additional physical constraints on the fit parameters. 
Gamma-ray spectra are particularly troublesome because 
they may contain significant contributions from inverse 
Compton emission, neutral-pion decay and nonthermal 
bremsstrahlung. By introducing additional physical con- 
straints, one can usefully reduce the set of linear com- 
binations of these models that fit the gamma-ray data. 
In §9, we discuss some of these constraints and describe 
how they may be imposed. 

2. PARTICLE DISTRIBUTION FUNCTION 

The algorithms used to compute the photon emission 
impose relatively few limitations on distribution func- 
tions suitable for use in fitting. The spectral models de- 
scribed below are derived assuming a particle momentum 
distribution function that depends only on the magni- 
tude of the particle momentum and not its direction. For 
most practical applications, the momentum distribution 
function should depend on a reasonably small number of 
parameters and should be integrable by adaptive quadra- 
ture rules. 

In applications to date, we have used a nonthermal 
particle distribution function of the form 
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Fig. 1. — Sample electron momentum distribution function. At 
low energies, a Maxwellian thermal distribution (feT = IkeV) is 
shown. At high energies, two nonthermal distributions are shown, 
each with a high-energy exponential cutoff i? C utoff = 10 TeV. The 
dashed line shows a nonthermal distribution with V = 2 and with 
positive curvature (a = 0.05) above momentum p = 1 GeV/c. The 
other nonthermal component (solid line) has T = 2 and zero cur- 
vature (a = 0). 




Fig. 2. — Example of the effect of curvature on nonthermal 
spectra. The solid lines show synchrotron (S), inverse Compton 
(IC), neutral-pion decay (7r) and nonthermal bremsstrahlung (NB) 
spectra with T = 2 and no curvature (a = 0). The dashed lines 
show the same spectra with positive curvature (a = 0.05) above 
particle momentum p = 1 GeV/c. 



where 



/(P)= 0, P<E /c, 



(2) 



and where p = ^mv, and Eq = 1 GeV. The same 
functional form is used for both protons and electrons. 
The normalization parameter, A, represents the density 
of particles with momentum p — Eq/c and has units 
cm- 3 (GeV/c)- 1 . 

When a = 0, equation (1) describes a power-law dis- 
tribution in momentum with an exponential cutoff at 
pc w £?cutoff- When a ^ 0, the power-law exponent 
changes with momentum for p > Eq/c (see Figure 1). 
The effect is such that, for each factor of ten increase in 
the momentum above Eq/c, the spectral index changes 
by an amount equal to the curvature parameter, a. Posi- 
tive values of a cause the particle spectrum slope to flat- 
ten toward higher momenta, as shown in Figure 1. Posi- 
tive values are expected due to nonlinear behavior of the 
diffusive acceleration mechanism (see e.g. Bell 1987; El- 
lison & Reynolds 1991; Berezhko & Ellison 1999). The 
qualitative effect of a small positive curvature on the 
shape of the nonthermal photon spectrum is shown in 
Figure 2. 

Each of the photon emission models described below 
depends on parameters associated with the underlying 
momentum distribution of nonthermal particles. From 
equation (1), these parameters are T, a and -E C utoff ■ The 
normalization parameter for each photon emission model 
includes the normalization parameter, A, for the relevant 
nonthermal particle distribution function. 

The low end of the particle momentum distribution is 
dominated by a thermal Maxwellian as shown in Figure 
1. As described in §9, one can impose charge conser- 
vation in the injection mechanism to set the value of 
the cosmic-ray proton normalization, A p relative to the 
cosmic-ray electron normalization, A e . 

In computing nonthermal photon spectra, integrals 
over particle momenta generally include only relativis- 



tic particles with 7 3> 1; in practice, we use 7 > 10. 



3. SYNCHROTRON RADIATION 

From Blumenthal & Gould (1970), the total syn- 
chrotron power emitted per unit frequency from an en- 
ergetic electron (7 3> 1) spiraling in a magnetic held is 



sin a / v 

^emitted [v) = 5 F \ — 

m e c z 



(3) 



where e is the electron charge, m e is the electron mass, 
B is the magnetic held strength, a is the pitch angle 
between the electron's velocity vector v and the magnetic 
held vector B, v c is the critical frequency, defined as 
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and F(x), the first synchrotron function, is defined as 

nOC 

F(x) = x dSK B/3 (Z), (x>0), (5) 



where K 5 / 3 (£) is an irregular modified Bessel function. 
Note that equation (3) applies to frequencies well above 
the gyro- frequency, where the synchrotron spectrum may 
be regarded as continuous. 

Assuming a steady state electron distribution, we de- 
fine N(p, a) dp dQ as the density of nonthermal electrons 
with pitch angles a within a solid angle dO, and momenta 
p within dp. From Blumenthal & Gould (1970), 



dfi Pcmittod(^) N(p,a) 



(6) 



is the total synchrotron power received per unit volume 
and per unit frequency, integrated over pitch- angles. 

Combining equations (3), (4) and (6), and assuming an 
isotropic distribution of pitch- angles so that N(p, a) = 
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N(p)/An, we obtain 
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Expressing this result in terms of a photon emission rate 
per unit energy, equation (7) yields the differential emis- 
sivity spectrum in the form 

1 dp N(p) R . . ( 8 ) 



dwdt hm e c 2 LO J \wo7^ 
where w is the photon energy and the angular integral is 



R(x) = - C da sin 2 a F f — ) 
w 2 7 Vsina/ 



(9) 



Crusius & Schlickciscr (1986) showed that equation (9) 
may be expressed analytically in terms of Whittaker 
functions, 

R( X ) = _ [ W (x; 0, 4/3) Wfo 0, 1/3) 

- W(x; 1/2, 5/6)W(z; -1/2, 5/6)] , (10) 

where W(z; k, fx) = e- z ' 2 z^ +1 ' 2 U (z; ± + n - k, 1 + 2/j) 
and U(z; m, n) is a confluent hypergeometric function of 
the second kind. 

The synchrotron flux is obtained by evaluating equa- 
tion (8). To speed numerical computations, we com- 
pute R(x) using a cubic spline interpolation on a pre- 
computed table; we use the GNU Scientific Library 
(Galassi et al. 2005) to perform numerical integration 
and one-dimensional spline interpolation and to evaluate 
selected special functions. This table is constructed by 
evaluating equation (10) on an adaptive grid of x values 
chosen to accurately sample the behavior of the function 
over the range 10~ 38 < x < 100. A log-spaced grid of 
x values was refined by adding interpolation points until 
the interpolated value of logi?(x) had a fractional error 
smaller than 1.25 x 10~ 5 at the midpoint of each x in- 
terval; we linearly interpolated the logarithms of x and 
R(x) only while refining the grid for the lookup table. In 
the final table, the associated cubic spline interpolation 
errors are \5R/R\ <5x 10~ 10 over the entire range of x. 

The parameters of the synchrotron model are T, a, 
-^cutoff, the total magnetic field strength, B to t and the 
normalization, 



■A/s 



ird 2 J v 



dVA e (r) 



A e V s 



(11) 



4ird 2 J v ew Aird 2 ' 

where A e is the normalization of the nonthermal electron 
momentum distribution, d is the distance to the source 
and Vs is the synchrotron emitting volume. A homoge- 
neous emitting volume is assumed. 

Note that the parameters of the model describe the 
physical properties of the synchrotron-emitting plasma 
rather than the properties of the observed synchrotron 
emission. In some situations, two or more of the fit pa- 
rameters are degenerate. For example, when fitting X- 
ray observations alone, the data constrain the critical 
frequency, v Cl but, because v c depends on the product 
-Btot7 2 , the model parameters B to t and -E C utoff are degen- 
erate. Similarly, when fitting radio observations alone, 



the normalization, A/s, and Btot are degenerate. In both 
cases, freezing B to t solves the problem. Although one 
could group the physical parameters to remove degenera- 
cies for special cases, we have chosen to keep the physi- 
cal parameters separate. The advantage of this choice is 
that, by introducing additional observational constraints 
from other energy bands, it may be possible to constrain 
the physical parameters separately. For example, by fit- 
ting radio, X-ray and gamma-ray data simultaneously, 
one might individually constrain A/s, T, a, -E C utoff and 
Btot- 

4. INVERSE COMPTON SCATTERING 

Given a distribution of relativistic (7 S> 1) elec- 
trons N(p), immersed in an isotropic radiation field with 
photon number density n(w,), the differential emissiv- 
ity spectrum of Compton scattered photons for single- 
scattering is 

dn 
dujdt 



! / duJi n{uJi) J dp N{p)gkn{i, <*>), 

(12) 



where 10 = hv/(m e c r), and 

okn(7> <*>i, w ) = 7 
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& Gould 1970), where 
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2(1 + Tq) 

is the Klein-Nishina scattering cross-section (Blumcnthal 



(14) 



4^7(7 - LU) 



T = 4cji7 and r = e 2 / (m e c 2 ) is the classical electron 
radius. 

The allowed range for q follows from the kinematics of 
Compton scattering. In the frame in which the electron 
is initially at rest, the incident photon has energy and 
the scattered photon energy may span the range oj' i /(l + 
2(1^) < ui' < In the lab frame, the scattered photon 
energy lies in the range u>i < uj < T/(l + T). Therefore, 
the allowed range for q is 



< q < 1. 



(15) 



limn - — - , / 

47(7 - u t 

The lower limit of the integral over electron momenta 
in equation (12) corresponds to the minimum electron 
momentum, p m i n = 7 m i n m e w, that can Compton scatter 
a photon from initial energy Ui to final energy ui. From 
the kinematics, one can show that the threshold electron 
Lorentz factor is 



7n 



(16) 



which corresponds to q = 1. 

The inverse Compton model is obtained by evaluating 
equation (12) for the cosmic background radiation field 
so that 



n(u>i) = 



1 



r 2 A 3 g^/e _ 1 ' 



(17) 



where A = h/ (m e c) is the electron Compton wavelength, 
9 ee kT/(m e c 2 ), and T = 2.725K (Bennett et al. 2003). 
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For computational convenience, it is useful to change 
the order of integration. Performing the integral over 
photon energies first, we let the electron momentum in- 
tegral extend over the full range of electron momenta and 
use the cross-section 



<7(7,Wi,w) 



CRN (7, Wj, w), 

0, 



<7min < 3 < 1; 

otherwise. 



(18) 



We can then speed up numerical computations by tabu- 
lating the photon energy integral, 



7 c (w,7) = J duji n(wi) o- KN (^,uj i} uj) 



(19) 



for a given spectrum of seed photons, n(wj). In practice, 
we use q as the variable of integration to evaluate the 
integral in equation (19). 

Because the range of I c (uj,j) spans many orders of 
magnitude over the domain of interest, we find it use- 
ful to interpolate on the logarithm of the function value. 
In constructing an interpolation table for I c (u>, 7), we ex- 
plicitly incorporate the existence of the threshold Lorentz 
factor 7 m ; n . Because I c (uj 7 j) asymptotically goes to 
zero as 7 — > 7 m i n , we set 7 c (o>,7) = for 7 < 70 = 
7min(wj jmax ) where Wi imax is the maximum incident pho- 
ton energy of interest. Introducing a change of vari- 
ables, x c = log (log (7/70)) and y c = logw, we tabulate 
log7 c (x c , y c ) on a 1024 x 1024 uniform rectangular grid 
in x c and y c covering the range 10 < 7 < 10 9 ' 5 and 
10 2 < uj < 10 15 eV. The smoothness of the resulting ta- 
ble allows accurate two-dimensional interpolation with a 
6 th -order spline (de Boor 1978) in each coordinate. 

The parameters of the inverse Compton model are T, 
a, E cuto f[, the blackbody temperature, T and the nor- 
malization, 
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where V\ is the homogeneous emitting volume. By de- 
fault, the temperature parameter, T, is not used and 
the model uses a lookup table for I c (uj,j) appropriate 
for seed photons from the cosmic background radiation 
field. The error associated with interpolation in this ta- 
ble is typically \SI C /I C \ < lO" 10 . 

A switch is provided to force the model to compute 
the radiation field integral by direct integration instead of 
table interpolation. The computational expense of direct 
integration currently makes this mode impractical for use 
in spectral fitting. 

Alternatively, one can use a linear combination of 
lookup tables to describe a more complicated radiation 
field. For example, one can construct lookup tables corre- 
sponding to a dilute stellar radiation field or to emission 
from molecular clouds. During the fit, one can vary the 
relative proportions of these components. Note that the 
current implementation restricts the input radiation field 
shape to one that is integrable by adaptive quadrature 
rules. 

The inverse Compton emitting volume, V\, need not be 
the same as the synchrotron emitting volume, Vs. Be- 
cause the cosmic background radiation photons will fill 
the entire synchrotron emitting volume, V\ will usually be 
at least as large as Vs. However, if nonthermal electrons 
are found in a volume with a relatively weak magnetic 



field, that volume will produce inverse Compton emis- 
sion, but little synchrotron emission, and V\ > Vs. For 
this reason, it is useful to allow the inverse Compton and 
synchrotron norms to be different. 

5. NONTHERMAL BREMSSTRAHLUNG 

In this section, we consider nonthermal brcmsstrahlung 
emission from a population of nonthermal electrons in- 
cident on a stationary target containing free electrons 
and ions. The total bremsstrahlung emissivity is com- 
puted as a sum of contributions from electron-electron 
and electron-ion bremsstrahlung. 

Given two populations of relativistic particles with mo- 
mentum distributions, N\(pi) and N 2 (p2), and with in- 
teraction cross-section, a, the general expression for the 
collision rate per unit volume is 



dn 
~dt 



= (1 + <*12) 1 J dpiJVi(pi) J dp 2 N 2 (p 2 



x ^(vi-vjf-fvjxvjf/c 2 , (21) 

where the integrals extend over the momenta of the in- 
teracting distributions (Landau & Lifshitz 1975). The 
Kronecker delta, S12 corrects for double-counting when 
the particles are identical. 

When |v 2 | <C |vi|, equation (21) reduces to the famil- 
iar non-relativistic form, and one of the populations may 
be treated as a stationary target. For example, when the 
target particles may be characterized by a Maxwellian 
thermal distribution, the thermal motions of the target 
particles may be neglected as long as kT 2 <C (7 — l)mc 2 . 

5.1. Electron- Electron Bremsstrahlung 

In the limit that |v 2 | <C |vi|, equation (21) yields a dif- 
ferential emissivity spectrum for electron-electron brems- 
strahlung of the form 



dn 
dcodt 



dp N(p) v 



dcr ee 
duj 



(22) 



where n e is the target electron density, and do~ ee /duj is 
the lab-frame differential cross-section for this interac- 
tion between identical particles and includes the factor 

(i + s 12 y\ 

The lab- frame cross-section, do~ ee / 'dujdip, differential in 
both photon energy, uj, and photon emission angle, tp was 
taken from Haug (1975) and was computed using soft- 
ware kindly provided by E. Haug. The lab- frame cross- 
section, differential in photon energy, dcr ee /duj 7 is ob- 
tained by integrating over photon emission angles. The 
angular integration limits are given by Haug (1975) and 
follow from energy-momentum conservation. In the cen- 
ter of momentum (CM) frame, the photon emission is 
symmetric along the direction of motion and all values 
of ip* are accessible. For an electron incident with mo- 
mentum p = "/mv in the lab frame, beaming restricts 
photon emission angles to a narrow cone in the forward 
direction. From Haug (1975), photons with energies 
(7 - l)/( 7 + 7/3 + 1) < uj < (7 - l)/( 7 -7/3+I) are 
emitted into a cone with the maximum emission angle, 
V'max, given by 



COSf/'max = 



(7 + 1)^ - (7 - 1) 
ury/3 



(23) 
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Lower energy photons may span < ip < ir. 

In carrying out the angular integration for ultra- 
relativistic electrons, some care must be taken to mini- 
mize numerical problems due to round-off errors. In par- 
ticular, the expression for the differential cross-section, 
d<r ee /dujdtp, includes terms that divide by the quantity 



x = lu(E — pcosip), 



(24) 



where E = T + l,p= ^T(T + 2), and T is the incident 
electron kinetic energy in units of m e c 2 . For 7 > I0 8 , x is 
identically zero in double-precision for cos^ = 1, causing 
a division by zero error. In the ultra-relativistic regime, 
we computed x using the first few terms of its series ex- 
pansion in powers of 7 -2 . Substituting cos-0 = 1 — s, and 
using s as the variable of integration, the relevant cancel- 
lations can be handled analytically. To check the result, 
we used the Lorentz invariant, Lod 3 a/dp 3 , to transform 
Haug's CM-frame cross-section into the lab frame (see 
e.g. Dermer 1986). In the CM frame, the cross-section 
computation is less affected by round-off errors because 
the relevant electron Lorentz factor is j c = ((7+l)/2) 1 / 2 . 
Applying the Lorentz transformation and treating round- 
off errors by handling the relevant cancellations ana- 
lytically, we verified that the transformed result agreed 
with the lab-frame cross-section in the ultra-relativistic 
regime. 

To speed up the numerical integrations, the cross- 
section is evaluated using a two-dimensional cubic-spline 
interpolation on pro-computed tables. In constructing 
an interpolation table for da ee /dcj = a', we explicitly 
incorporate the fact that cr'(T, u>) = for kinetic ener- 
gies T < T m i n (u>), where T^^lo) is determined by the 
kinematics. We define T m j n (cj) numerically as the lo- 
cus of points at which a' (T m ; n , u>) ~ eBmax{er'(T, u>)} 
for = 1CT 8 . Introducing a change of variables, 
xb = log(log(T/T m ; n )) and ys = logw, we tabulate 
log <t'(xb, Vb) on a 1024 x 1024 rectangular grid in xb and 
ys covering 10 2 < T/mc 2 < 10 95 and 10 2 < uj/mc 2 < 
10 9 5 . The smoothness of the resulting table allows ac- 
curate two-dimensional cubic spline interpolation. The 
error associated with interpolation in the cross-section 
table is \Sa'/a'\ < 10" 5 . 

5.2. Electron-Ion Bremsstrahlung 

From equation (21) the differential emissivity spectrum 
for electron-ion bremsstrahlung can be written 



dn 
dujdt 



nz / dp N(p) v 



da 



cZ 



dui 



(25) 



where nz is the density of target ions with charge Z, and 
da e z/duj is the lab- frame differential cross-section. 

The Bethc-Hcitlcr cross-section (Hcitlcr 1953; Koch & 
Motz 1959) determines the probability that deflection of 
a relativistic electron in the unscreened field of an ion of 
charge Z will yield a photon of energy u> = hv/{m e c 2 ). 
The differential cross-section may be written in the form 



dcr e z _ 4> 7/3 
duj uj 70 (3 



+ 
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where 



x = 



^ 2 (1 + /OT , 



3/3o/3 7o7/3g/3 3 2j j(3o(3 
«o(7 + 7o/?o) a (7 o +lP 2 ) ^ _ 



7o7/3o 2 /? 2 



7o/?o 7 2 /? : - 70/Mo 

(p = aZ 2 r 2 , r = e 2 / (m e c 2 ), 



. (27) 
(28) 



a = 21n[ 7o (l+A>)], a=21n[ 7 (l + /3)], (29) 



L = 2 In [(707 + 7o7/?o/3 - 1) /w] . 



(30) 



In these expressions, a is the fine structure constant and 
7(/3) and 70(A)) are the Lorentz factors of the scattered 
and incident electron, respectively. The Bethe-Hcitler 
cross-section is derived using the Born approximation 
which is appropriate in the limit of high kinetic ener- 
gies, such that 2naZ/ ' [3 <C 1 for both the incident and 
scattered electron. The recoil of the nucleus is neglected 
so 7o(/?o) = j(fi) The accuracy at low energies is im- 
proved by including the Elwert correction factor, which 
we apply at all energies (Elwert 1939; Pratt & Tseng 
1975; Haug 1997), so that da e z/duj — ► rjEdcr e z / dui where 



and where 



_ £ 1 - exp(-Co) 

r l E = T~\ 7 — pr> 

Co 1 - exp(-£) 



Z Z 
£ = 2ira-, £, = 2ira— . 

P Po 



(31) 



(32) 



When the target ions retain bound atomic electrons, it 
is necessary to modify the cross-section to account for 
screening of the nuclear charge. We assume that the 
target material is completely ionized so that screening 
corrections are not necessary. 

The relative simplicity of the cross-section makes it 
practical to perform the integration over electron mo- 
menta by evaluating (26) directly rather than interpo- 
lating values from a pre-computed table. 

The parameters of the nonthermal bremsstrahlung 
model are T, a, -E C utoff and the normalization, 



^ J v dV n (r)A e (r) = 



(33) 



where no = ^2 z Uz 1S the total ion number density 
of the target, nz is the number density of ions with 
charge Z and Vb is the emitting volume. The user in- 
terface includes parameters to control the relative contri- 
butions of electron-electron and electron-proton brems- 
strahlung. The electron-electron contribution has weight 
Xe = J2z Znz/n and the electron-proton contribu- 
tion has weight Xi = ^ z Z 2 nz /n^. By default, the 
target is assumed to consist of hydrogen and helium 
with nHc/^H = 0.1 so that the default weights are 
X e = X H + 2A Hc = 1.091 and X l = X H + 4A Hc = 1-273. 
A common alternative is to view n as the target pro- 
ton density and to compute the weights using that as- 
sumption. As long as the weights are consistent with the 
definition of no, the result is the same. 
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6. NEUTRAL PION DECAY 

Nonthermal protons produce gamma-ray emission pri- 
marily through collisions with thermal protons. These 
collisions yield neutral pions via pp — > 7r° + X, and 
the neutral pions decay via ir° — > 27. We include only 
the contribution from proton-proton collisions and ignore 
contributions from processes other than the decay of neu- 
tral pions. In its rest-frame, the neutral pion decays 
within ~ 10 -16 sec, producing two gamma-rays, each 



with an energy of luq = 



Because the pion has 



zero spin, the gamma-rays are emitted isotropically. Fol- 
lowing Hillier (1984), the number of gamma-rays emitted 
into a rest frame angle between (9* and + d9* is 



n(0»)d0* = sin6>*d#» 



(34) 



Since a gamma-ray that is emitted at an angle 0* in the 
pion rest frame has a lab-frame energy 



UJ = Wo7tt (1 — Pir cos 0*) . 



(35) 



the number of gamma-rays with lab-frame energy be- 
tween lj and lj + du> is 



n{uj) = n(0») 



1 



d0. 



dw wo77r/5 



(36) 



where p v = 7 7I .m 7r w 7r is the pion momentum. For a pop- 
ulation of pions with momentum distribution, N- n {p 7r ), 
the gamma ray spectrum is then 



n{u) = 2 



Ptt.t, 



AT I \ d P* 
Pit 



(37) 



The lower integration limit is set by the minimum pion 
momentum required to yield a lab-frame photon of en- 
ergy lj. From equation (35), the low energy gamma-ray 
has energy w m i„ = wo7tt(1 — At), corresponding to a pion 
velocity of 



2 1 

P"7r,min — 9 , 9 • 



(38) 



A lab frame photon of energy lo therefore requires a pion 
Lorentz factor of at least 



Ttt,: 



1 / LO 

2 I 17 



LOq 

LJ 



(39) 



The upper integration limit corresponds to the maximum 
pion momentum that can be produced by a proton of 
kinetic energy T p ^ max . From the collision kinematics, one 
can show that 



1 (T, 



Ttt,: 



p.max 

2lj 



TO B C 2 



(40) 



The production spectrum of secondary pions in proton- 
proton collisions is 

NM = n p T dpv p N p (p) d<7{ ^' p) , (41) 

where n p is the density of target protons, v p is the non- 
thermal proton velocity, N p (p) is the nonthermal proton 
momentum distribution, and dcr{p n ,p)/dp„ is the dif- 
ferential cross-section for production of a neutral pion 
with lab-frame momentum p v from a proton with lab- 
frame momentum p. The lower integration limit is the 
threshold proton momentum for producing a pion with 




Fig. 3. — Distribution of fractional recurrence relation errors. 
The fraction, /, of errors larger than p, is shown for each spectral 
model. The errors were computed at many points along each of 
several model spectra. One spectrum was generated for each pair 
of values, T and E C utoffi over a grid spanning 1.8 < T < 4 and 
1 < ^cutoff < 1000 TeV. The dotted line corresponds to /(> p) = 
exp(— p 2 /2<T 2 ), with a = 2.5 X 10~ 5 . Note that the curves for 
electron-electron and electron-proton bremsstrahlung are almost 
identical. 



momentum p v . From the collision kinematics, one can 
show that the threshold proton kinetic energy is 



T, 



I \ o /„2 2 i 2 4\ V 2 . m -K C 
,, n.m(P7r) = 2 {p~(? + m%e) +"^- 



(42) 



For the purpose of interactive spectral fitting, com- 
puting the photon spectrum from neutral-pion decay by 
direct integration over the proton momentum distribu- 
tion and over the resulting pion spectrum is quite com- 
putationally demanding. The cost of the numerical inte- 
grations is compounded by the fact that the differential 
cross-sections are computationally expensive for certain 
energies (see e.g. Dermer 1986; Mori 1997). To reduce 
the cost of these computations, we evaluate the integral 
over the proton distribution using the delta-function ap- 
proximation described by Aharonian & Atoyan (2000). 

The parameters of the neutral-pion decay model are T, 



a, E r 



itofr 



and the normalization. 



i 



47rd 2 



dV n p (r)A p (r) 



And 2 ' 



(43) 



where A p is the normalization for the nonthermal pro- 
ton momentum distribution, n p is the density of target 
protons, and is the emitting volume. Note that for 
this model, the parameters T, a and -E C utoff refer to the 
proton momentum distribution. 



7. ACCURACY OF COMPUTED SPECTRA 
Consider spectral models of the form 
-r 

S r (u:;e 



dp(^ 



ex Pl f°— E ) <f,(p, u ), (44) 



where lo is the photon energy, p is the particle momen- 
tum, eg is a scaling constant, and cj>(p, lu) depends on the 
physical process. Taking the derivative with respect to 
l/e c = a, we find that 



1 dSr(uj;a) 
£o da 



Sr(w,a) - 5r-i(w;a). (45) 
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This identity specifies a recurrence relation between spec- 
tra of the form shown in equation (44) . 

Numerical evaluation of equation (44) yields model 
spectra of the form Sr(w;a) + A.(u>;a, T), where Sr is 
the exact result and A represents the error in the compu- 
tation. Because accurately computed spectra must sat- 
isfy equation (45), it follows that A(ui;a,T) must also 
satisfy equation (45). Although an error in the compu- 
tation of <f>(p,ui) might generate such (a,T) dependent 
errors in the computed spectrum, such an error should 
be detectable on comparison with an analytic solution 
or with an independent numerical calculation. Eliminat- 
ing that possibility, the error term in computed spec- 
tra that satisfy equation (45) must be independent of 
a and T so that A(oj;a,T) = A(u>). The magnitude of 
the energy-dependent error, A(w) must be estimated by 
other means, such as by comparing with analytic solu- 
tions and with independent numerical results. Although 
satisfying the recurrence relation does not prove that the 
spectral computations are correct, it does provide strong 
constraints on the magnitude and parameter dependence 
of the error term, A. 

To verify equation (45) numerically, we introduce a 
finite-difference approximation for the derivative. It 
follows that computed spectra of the form shown in 
equation (44) should obey a relationship of the form 
1Z(qj; r, a) = where 



1.02 - 



K{uj;T,a) = 1 



T-i 



a) 



— lim 

Sa^O 



S r (a) 
S r (a + 5a/2) 



S r (a-5a/2) 



Sr(a)eo5a 



(46) 



For clarity in equation (46), we have avoided writing 
out the explicit dependence of the right-hand side upon 
photon energy, to. The centered difference used to ap- 
proximate the derivative in equation (46) should yield 
quadratic convergence in the limit that 5a — > 0. For 
each of our spectral models, we verified that, in the limit 
5a — ► 0, 1Z — ► as 1Z cx (5a) 2 ; observing smooth conver- 
gence at the expected rate confirms that the models be- 
have as expected. In subsequent evaluations of equation 
(46), we adopt a value of 5a/ a = 2.5 x 10~ 6 . Because in- 
dividual terms in equation (46) can be ^ 1, it is useful to 
compute the fractional error, p = 1Z/\1Z\, where the de- 
nominator is the £ 2 -norm of the three terms in equation 
(46). 

To test the T and a dependence of A, we computed 
values of p at many points along each of several spec- 
tra. One spectrum was computed for each pair of (T, a) 
values on a grid sampling the range 1.8 < T < 4 and 
1 < -Ecutoff < 1000 TeV. For the synchrotron model, p 
was evaluated at photon energies in the range 10 -5 eV < 
uj < 10 5 w c , where uj c is the critical energy. For the other 
models, p was evaluated at photon energies in the range 
10 5 eV < uj < 3£ , cu toff- Each photon energy grid was log- 
arithmically spaced and used 80 points per decade. Fig- 
ure 3 shows that, throughout this range, p was < 10~ 4 . 
We conclude that, throughout the parameter range of in- 
terest, the (a, T) dependence of A is no larger than about 
a part in 10 4 . This result suggests that the errors are es- 
sentially independent of the parameters of the particle 
distribution function. 

To test the dependence of the errors on the other pa- 
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Fig. 4. — Ratio showing our synchrotron models divided by those 
of Sturner et al. (1997) for T = 1.8 (dashed), T = 2 (solid), and 
r = 2.3 (dotted). Note that the three curves are almost identical. 



rameters, including the photon energy, we compare our 
results with other analytic and numerical solutions. Blu- 
menthal & Gould (1970) derive the well-known analytic 
result that the synchrotron spectrum from a power-law 
distribution of electrons, N^) cx 7~ r , is itself a power- 
law of the form S(y) cx i)/ 2 . To compare with 
this result, we used a power- law electron distribution 
function to compute numerically S , syn c(^) = Lodn/dcudt, 
as described above in equation (8), and we computed 
S(v) = dW/dvdt from equation (4.59) of Blumenthal 
& Gould (1970). In deriving the analytic solution, the 
lower limit of the integral over 7 is extended to zero (Blu- 
menthal & Gould 1970), effectively neglecting the elec- 
tron mass. Consistent with this assumption, the value 
of x = vjv c in the numerical integrand must be com- 
puted using the approximation 7 « p/(m e c) otherwise, 
the computed spectrum departs from the analytic solu- 
tion at low frequencies. We find that the absolute value 
of the fractional error is 



^sync I 



1 



S{y) 



< 3 x 10 



-11 



(47) 



for frequencies in the range 10 7 < v < 10 20 Hz and 
for values of T and i?tot in the ranges 1 < T < 4 and 
1 < -B tot < 10 4 pG, respectively. The primary source 
of error in this comparison is associated with interpola- 
tion of R(x) in our precomputed table; computing R(x) 
directly in terms of hypergeometric functions, we repro- 
duce the analytic solution to within |e syn c| < 3 x 10~ 13 
over the specified range of v, T and -B to t- 

Blumenthal & Gould (1970) also discuss the inverse 
Compton spectrum produced by a power-law distribution 
of electrons scattering photons from a blackbody radia- 
tion field. They give asymptotic analytic solutions valid 
in the Thomson limit and in the extreme Klein-Nishina 
limit. To compare our inverse Compton model with these 
analytic solutions, we used a power-law electron distri- 
bution function, ^(7) cx 7~ r , to compute numerically 
Sinvc(w) = dn/dtodt, as described above in equation (12). 
For the purpose of these tests, we evaluated the integral 
over incident photon energies (equation 19) by direct nu- 
merical integration rather than by spline interpolation in 
precomputed tables. 

In the Thomson limit, the energy of the incident pho- 
ton in the electron rest frame is small compared to the 
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Fig. 5. — Ratio showing our inverse Compton and nonthermal 
bremsstrahlung spectra divided by those of Sturner et al. (1997) 
for T = 1.8 (dashed), T = 2 (solid), and T = 2.3 (dotted). The 
three curves for the inverse Compton spectra are almost identical. 



electron rest energy and the momentum transferred by 
the electron is small compared to its initial momentum. 
In this regime, corresponding to < w < jmc 2 , the 
scattering cross-section is essentially independent of the 
energy of the incoming photon. The Thomson limit ap- 
proximation is most accurate in the limit of low radia- 
tion temperature, T, and for scattered photon energies 
near the middle of the applicable energy range. Such a 
case provides the best test of the accuracy of our inverse 
Compton model in the Thomson limit. For power-law 
slopes in the range 1 < T < 4 and radiation field temper- 
atures 10 -2 < T < 10 4 K, we used our model to compute 
inverse Compton spectra, Si nvc (o;), for scattered photon 
energies in the range 107^ in /cT < u < 0.l7 m i n mc 2 , us- 
ing 7 m ; n = 10. For comparison, we computed s(lj) = 
diVtot/didei oc uj- {t+1 ^ 2 using equation (2.65) of Blu- 
menthal & Gould (1970). We found that our inverse 
Compton model converged smoothly to the analytic re- 
sult in the appropriate limit. The smallest fractional 
difference is 



1 



< lo- 



rn 



Broader comparisons are less useful as a test of com- 
putational accuracy because the asymptotic solution it- 
self becomes less accurate with increasing T and for 
energies approaching the endpoints of the applicable 
range. For T = 0.01 K, the fractional error is smallest 
(|einvc| < 10 -8 ) at a; 50 eV, then increases toward lower 
and higher energies with |ei nV c| < 10 -3 for the range 
nr 2 < lo < 10 7 eV. Increasing T narrows the applicable 
energy range so that for T — 10 4 K, |ei nvc | < 10 -3 only 
for the range 2 < uo < 5 keV. 

In the extreme Klein-Nishina limit, the energy of the 
incident photon in the electron rest frame is large com- 
pared to the electron rest energy. In this regime, cor- 
responding to LOi^i 3> mc 2 , the scattering cross-section 
is strongly peaked near the maximum scattered pho- 
ton energy so that individual Compton scatterings tend 
to involve a large energy transfer. The characteris- 
tic scattered photon energy is then ui ~ -fmc 2 . It 
follows that this regime may be characterized by the 
requirement that uiiuj 3> m 2 c 4 so that, for a black- 
body radiation field with Ui ~ kT, the extreme Klein- 



Nishina limit corresponds to scattered photon energies 
u > m 2 c 4 /(fcT). We computed s(u>) = djV tot /dtdei us- 
ing equation (2.88) of Blumenthal & Gould (1970) 1 . In 
the extreme Klein-Nishina limit, we computed inverse 
Compton spectra, 8i nvc (w), for scattered photon energies 
in the range 10 3 m 2 c 4 /(fcT) < u < 10 12 m 2 c 4 /(/fcT) for 
power-law slopes in the range 1 < V < 4 and radia- 
tion field temperatures in the range T = 10 3 — 10 8 K. 
We found that the inverse Compton model spectra con- 
verged smoothly to the analytic result, s(u>). The frac- 
tional error decreased smoothly with increasing photon 
energy from |e in vc| ~ 10" 3 at w « 10 3 m 2 c 4 /(fcT) to 
|£invc| ~ 3 x 10~ 8 for u « 10 9 m 2 c 4 /(fcT). Numerical 
round-off errors become important at higher energies. 

To examine the absolute accuracy of our computed 
photon spectra for more realistic particle spectra, we 
compared our results with those of Sturner et al. (1997). 
To facilitate this comparison, Sturner kindly computed 
photon spectra for a particle distribution function of the 
form 



N(p) = A' ( 



pc 



1 McV 



exp 



T 



toir 



(49) 



where A' has units of cm~ 3 (MeV/c)~ 1 and T is the 
kinetic energy. Sturner provided synchrotron, inverse 
Compton and nonthermal bremsstrahlung spectra for 
T = 1.8, 2.0 and 2.3 and E cutoS =10 TeV. The syn- 
chrotron spectrum was computed using B tot =l /iG. The 
inverse Compton spectrum was computed for the cos- 
mic background radiation field, using a temperature of 
2.7 K. The nonthermal bremsstrahlung spectrum was 
computed for a fully ionized target with ion den- 
sity 0.11 cm -3 consisting of protons (0.1 cm~ 3 ), alpha- 
particles (0.01cm -3 ), and free electrons (0.12cm -3 ), 
corresponding to relative weights of X e = 1.090 and 
Xi = 1.182, for electron-electron and electron-ion brems- 
strahlung, respectively. We used our models to compute 
photon spectra for the same parameters and particle dis- 
tribution function. 

Figures 4 and 5 show our spectra divided by those ob- 
tained from Sturner. Over most of the energy range, 
our spectra agree quite well with Stumer's; the inverse 
Compton and synchrotron spectra agree to within < 1%, 
while the nonthermal bremsstrahlung spectra typically 
differ by < 5%. Aside from the weak T dependence in 
the nonthermal bremsstrahlung differences below about 
100 MeV, the differences between Stumer's spectra and 
ours are largely independent of the spectral index, T. 
Because the T dependent errors in our spectra are con- 
strained by the recurrence relation, the weak residual T 
dependence seen in the nonthermal bremsstrahlung dif- 
ferences appears to be associated with Stumer's spec- 
tra. The reason for the overall ~6% discrepancy be- 
low 0.25 MeV is unclear; features in the nonthermal 
bremsstrahlung ratio near 30 TeV and near 10 MeV cor- 
respond to points in Stumer's spectra at which the slope 



1 Equation (2.88) of Blumenthal & Gould (1970) involves a con 
stant, C;, that is defined in terms of an infinite series that con- 
verges extremely slowly: C; 



(6/tt 2 ) Efcilnfc/fc 2 ~ 0.5700. 
From its series definition, it follows that this constant is equiva- 
lent to C; = — (6/71- 2 ) C'(2), where f is the first derivative of the 
Riemann zeta function. To enable more precise quantative com- 
parison with the analytic solution, we used the numerical value 
Ci = 0.56996099309453280, obtained using the symbolic algebra 
package MAPLE. 
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Fig. 6. — Ratio showing our synchrotron fluxes divided by SRCUT 
synchrotron fluxes from XSPEC for T = 1.8 (dotted), T = 2 (solid), 
T = 2.3 (dashed). We used B tot = 10 (J.G and E cutoB =10TeV 
corresponding to inbreak — 
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Fig. 7. — Ratio showing our neutral-pion decay spectrum divided 
by that of Mori (1997). 



changes discontinuously. Because Stumer's synchrotron 
and inverse Compton fluxes become negative at the ex- 
treme end of the high-energy photon spectra, rather than 
asymptotically approaching zero as ours do, we attribute 
the larger differences near the cutoff in these spectra to 
numerical errors in Stumer's spectra. However, since 
the largest differences occur at extremely low flux lev- 
els, they are probably not important for observational 
comparisons. 

We also compared our synchrotron model with the SR- 
CUT model (Reynolds 1998; Reynolds & Keohane 1999) 
in xspec version 11. 3. 2d. SRCUT is designed to com- 
pute the synchrotron spectrum from an exponentially 
cut off power-law distribution of electrons in a homo- 
geneous magnetic field. It depends on a break-frequency 
parameter, ^brcak, defined as the critical frequency for 
sin a = 1 and 7 = E cuto g /(m e c 2 ). For v < ^break, we 
verified that SRCUT is consistent with the analytic result 
for a power-law electron distribution. For frequencies 
comparable to ^brcak or larger, we find that SRCUT spec- 
tra differ from our spectra by as much as a factor of two 
or more. Figure 6 shows that the ratio of the two models 
depends on photon energy and on the power-law index. 
The ratio also depends on the break frequency in the 
sense that the models agree for v <C i-Wak but disagree 
for frequencies near the break and above. We are con- 
fident that our computed synchrotron spectrum is accu- 
rate, first, because the validity of equation (45) indicates 
that any errors are independent of the power-law index, 
r, and cutoff-energy, i? C utoff and, second, because our 
computed spectrum is consistent with Stumer's model. 
We conclude that, for frequencies near and above inbreak, 
the SRCUT model does not accurately represent the syn- 
chrotron spectrum from an exponentially cut-off power- 
law distribution of electrons in a homogeneous magnetic 
field. 

We tested our neutral-pion decay model by comparing 
our results with the results of Mori (1997). For the pur- 
pose of this comparison, we used his proton momentum 
distribution function of the form N(p) — (4t:/v) J p (p) 



where 



J P (p) = ' 



6.65 x 10" 



-2.75 



.100 Gcvy 



E > 100 GeV; 



1.67 



(-*- 

\GcV/<l 



-2.7 



^ 2.5GeV/ c y 



-1/2 



where p = jrripV and E = T + m p c 



E < 100 GeV, 

(50) 

2 . Figure 7 shows our 
computed neutral-pion decay gamma-ray spectrum di- 
vided by the fluxes given in Table 1 of Mori (1997). The 
largest difference occurs near 100 MeV, where our flux is 
about a factor of two larger than that of Mori. Because 
Mori (1997) used a much more detailed model of pion 
production, we do not expect exact agreement. Yet, for 
photon energies above 1 GeV, the two gamma-ray fluxes 
agree to within about 10%. Because we are primarily 
interested in fitting data in the 1-10 TeV band, our sim- 
plified pion-decay model is adequate for our needs. 

8. FITTING OBSERVED SPECTRA 

Fitting a model to an observed spectrum involves min- 
imizing a goodness-of-fit statistic that compares the ob- 
served spectral data to the spectral model. This com- 
parison often involves binned data. 

For example, consider X-ray observations which pro- 
vide the observed number of counts in each energy bin. 
Neglecting nonlinear effects in the X-ray detector, the 
expected number of counts is usually computed using an 
expression of the form (Davis 2001) 

C{h)=B{h)+t [ dE R(h,E)A(E)S(E), (51) 

where C(h) is the total number of counts in bin h, B(h) 
is the number of counts due to the instrumental back- 
ground, t is the exposure time and S{E) is the spectral 
model describing the incident flux of photons with en- 
ergy E. In equation (51), R(h,E) is the redistribution 
function, describing the probability that incident pho- 
tons with energy E contribute counts to bin h, and A(E) 
is the effective area, accounting for the telescope collect- 
ing area, the transmission efficiency of the optical system 
and the quantum efficiency of the detector. 
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In spectral fitting software, equation (51) is usually 
implemented as a discrete sum of the form 

C h = B h + tY J Rh,kA k S k . (52) 

k 

The redistribution function, R(h,E), is represented as a 
matrix, Rh,k, and the effective area and source models 
are represented as vectors, Ak, and Sk, respectively. 

To accurately represent the integral in equation (51), 
the software must compute the model, Sk, as an integral 

S k (E k ;AE k ) = f dES(E), (53) 

JAE k 

over the width AE k of each spectral bin, E k . For binned 
data, we evaluate these integrals using the well-known 
Simpson's rule. Although this approach requires three 
function evaluations per bin, the bin-edge function values 
may be shared between neighboring bins, reducing the 
total cost of each spectrum computation by about 30%. 
This approach is accurate to the extent that a quadratic 
polynomial is a good approximation to the underlying 
function, S(E), within each spectral bin. 

9. FIT CONSTRAINTS 

When simultaneously fitting radio, X-ray and gamma- 
ray observations of supernova remnants, the degeneracy 
of certain fit parameters (§3) and the variety of emission 
mechanisms in the gamma-ray band may make it diffi- 
cult to determine a unique set of fit parameters unless 
additional constraints are available. 

Based on reasonably general principles, a number of 
constraints can sometimes be imposed. For example: 

1. When the magnetic field associated with the syn- 
chrotron emission is generated primarily by cosmic- 
ray streaming (Lucek & Bell 2000), the energy den- 
sity in cosmic-rays should set an upper limit on the 
magnetic energy density. The corresponding upper 
limit on the magnetic field strength is 

Stot < y^Tr (e e +e p ), (54) 

where e e and e p are the energy density in cosmic- 
ray electrons and protons respectively. Because 
£?tot and -E cu toff are degenerate, introducing an up- 
per limit on B to t effectively sets a lower limit on 

-Ecutoff • 

2. If electrons and protons are injected into the accel- 
erator at the same rate, the normalization of the 
proton momentum distribution A p can be fixed by 
requiring equal densities in nonthermal electrons 
and protons at some characteristic injection kinetic 
energy T inj < m e c 2 : 

N e {p e s n j)dp e = N p (p pM - 3 )dp p . (55) 

When both distribution functions are of the form 
N(p) = A(p/p )~ r and share the same power-law 
exponent, T, this constraint implies that A p /A e = 
(mp/TOe)^- 1 )/ 2 ss 91 if r = 2.2 (Bell 1978). By 
fixing the ratio A p /A e , this constraint reduces the 
variation in the ratio A/^/Ab to the variation asso- 
ciated with the mass ratio of the associated targets, 
n p V^/{n VB)- 



3. In many cases it should be reasonable to as- 
sume that the emitting volume that produces syn- 
chrotron emission will also produce inverse Comp- 
ton emission due to up-scattering of cosmic back- 
ground photons. In such cases, one can impose 
a lower limit on the normalization associated with 
this inverse Compton process so that A/^cbr > A/"s. 
This constraint ensures that fits to the gamma-ray 
spectrum will include an appropriate contribution 
of inverse Compton emission. A similar argument 
can be used to constrain the minimum value of the 
nonthermal bremsstrahlung normalization, Ab for 
an assumed minimum target density. 

To support imposing constraints based on charge con- 
servation, our software provides a function to compute 
the proton norm, A p , that, for a given electron norm, A e: 
will yield equal nonthermal electron and proton densities 
at a given injection kinetic energy Ti n j. For constraints 
that depend on the energy density in nonthermal parti- 
cles, our software provides functions that compute the 
energy density of each particle population. The energy 
density is defined to be 

/>oc 

e = mc 2 I dp N(p) ( 7 (p) - 1) , (56) 

J Pmin 

where the lower integration limit, p m i n , is the momentum 
at which the thermal and nonthermal particle densities 
are equal. Note that p m ; n depends on the density and 
temperature of the thermal particles and on the non- 
thermal particle momentum distribution. 

In practice, complicated fit constraints such as the up- 
per limit on B to t may be imposed using a technique anal- 
ogous to the method of Lagrange multipliers (Mathews & 
Walker 1965). Rather than minimizing % 2 , the idea is to 
construct a constraint function g(x) > which may de- 
pend on a vector of parameters, x, and then to minimize 
the sum, x 2 + -\9( x )> where A is a parameter that deter- 
mines the importance of the constraint. The constraint 
function should be constructed to ensure that g(x) = 
when the constraint is satisfied. For example, to impose 
the constraint that _B 2 ot < 87r(e e + e p ), one might choose 

aU) = { B tot [8^(£e + ep)}- 1 , Bl t > 87r(e e + e p ); 
yy ' \0, otherwise. 

(57) 

Similar terms may be added to impose additional con- 
straints. 

10. CONCLUSIONS 

We have described models that can be used to compute 
the synchrotron, inverse Compton, nonthermal brems- 
strahlung and neutral-pion decay spectra of homoge- 
neous sources containing nonthermal electrons and pro- 
tons. We have implemented these models as an im- 
portable module for use in the spectral fitting package 
ISIS (Houck & Denicola 2000); the software is available 
from the ISIS web page 2 . The models are designed to 
be accurate and fast enough for use in interactive data 
analysis on a typical workstation. To the best of our 
knowledge, this is the first implementation of some of 
these models in a form suitable for interactive fitting with 

2 http://space.mit.edu/cxc/isis 
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publicly available data analysis software. Results derived 
using the synchrotron and inverse Compton models have 
been presented elsewhere by Pannuti et al. (2003) and 
by Allen, Houck & Sturncr (2005). 

We assessed the accuracy of these models by using a 
recurrence relation, by comparing with analytic solutions 
for synchrotron and inverse Compton scattering and by 
comparing with published work by other authors. The 
accuracy with which our solutions obey the recurrence 
relation in equation (46) demonstrates their correct de- 
pendence upon the power-law index and cutoff-energy 
parameters. We also showed that, over most of the en- 
ergy range of interest, our models agree with those of 
Sturner et al. (1997). The largest differences between 
our models and those of Sturner et al. (1997), of order 
~10%, occur near the cutoff of each model spectrum and 
are essentially independent of the spectral index. 

We found much larger differences between our syn- 
chrotron spectrum and the SRCUT model from xspec. 
Although consistent with the analytic solution for v <C 
^broak, SRCUT differs from our model and from Stumer's 
by as much as a factor of two or more in the X-ray band 
near ^Wak (see Figure 6). The most important difference 
is that the normalization derived using SRCUT overesti- 
mates the radio flux at 1 GHz by an amount dependent 
upon the spectral index as shown in Figure 6. The spec- 

3 http://glast.gsfc.nasa.gov 



tral index from SRCUT tends to be a few percent too steep 
and the break frequency tends to be a few percent too low 
but, in practice, such differences may be detectable only 
with very high quality data. Note that, by overestimat- 
ing the radio flux, SRCUT may suggest the existence of 
positive curvature in the underlying particle momentum 
distribution. 

In future work, we hope to improve the neutral-pion 
decay model by explicitly computing the integral over 
proton momenta using improved pion-production cross- 
sections. This refinement will extend the applicable 
range to sub-GeV photon energies and will improve the 
overall accuracy of the model. These improvements will 
be more important as better observations of the gamma- 
ray spectrum become available from the Gamma-Ray 
Large Area Space Telescope (GLAST) 3 and from future 
advances in instrumentation. 



We thank E. Haug for providing software to compute 
values for the electron-electron bremsstrahlung cross- 
sections and S. Sturner for providing numerical tables 
used to verify our results for selected electron momentum 
distributions. We also acknowledge useful conversations 
with John E. Davis and a number of comments from the 
referee which helped us to improve the paper. 



APPENDIX 
A. SPECTRUM TABLES 

For reference, Tables Bl, B2 and B3 contain sample spectra for each emission process described in this paper. The 
nonthermal electron and proton momentum spectra have the same shape, with slope r = 2, and a cutoff energy 
^cutoff = lOTeV. Table Bl gives synchrotron spectra for curvatures a = and a = 0.05. Tables B2 and B3 give 
the inverse Compton, nonthermal bremsstrahlung and neutral-pion decay spectra for curvatures a = and a = 0.05, 
respectively. All normalization coefficients were set to unity. The synchrotron spectrum was computed for B to t = 
10 fiG. The inverse Compton spectrum was computed using a 2.725 K blackbody distribution of seed photons (Bennett 
et al. 2003). Contributions to nonthermal bremsstrahlung due to the electron-electron and electron-proton processes 
are listed separately, with their respective weights set to unity. 

B. DISTRIBUTED COMPUTATION 

The spectral models described in this paper were designed to fit multi-wavelength spectral data interactively on a 
typical workstation and to achieve a high degree of accuracy as efficiently as possible. In practice, to be sure that the 
fit has fully converged and to derive confidence limits, it is necessary to thoroughly examine the parameter space in 
the neighborhood of the best-fit parameters. To conduct this search more quickly, we have found it useful to distribute 
the task of computing single-parameter confidence limits over a number processors running in parallel (Noble et al. 
2005). One master process manages the computations being performed by a number of slave processes that run on 
different computers. All of the computers are connected by a local network and all processes have access to the relevant 
data and spectral models. The master process assigns each slave process the task of computing confidence limits for 
a single parameter. If a slave process finds an improved fit, that new parameter set is sent to the master process. If 
that fit is the best yet found by any slave, the master commands all the slave processes to re-start the search from 
the beginning, using the new parameter set. We have found that the time required to obtain a set of converged single 
parameter confidence limits using this approach is often reduced by more than a factor of N } where N is the number 
of slave processes. 

We have implemented this algorithm using the Parallel Virtual Machine (pvm) (Geist et al. 1994) to handle message 
passing between the master and slave processes. We used the spectral fitting package, ISIS (Houck & Denicola 2000) to 
perform model fits and confidence limit searches using a set of S-Lang scripts. An S-Lang module provides a scriptable 
interface to the PVM library, making it possible for these scripts to communicate with the pvm. 
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TABLE Bl 

Sample Synchrotron Spectra 



Flux Flux 





Energy 




a=0 




a=0.05 




Energy 




a=0 




a=0.05 




(eV) 




(photons s 1 


cm 


- 2 GeV- 1 ) 




(eV) 


(photons s 


cm 


" 2 GeV" 1 


1 


.OOOOe- 


-07 


6 


R9Q4p+1 R 


6. 


,8613e+15 


1 


,7783e-01 


2 




7. 


,0374e+06 


1 


,7783e- 


-07 


2 


. O 1 sOc T J. \J 


2 


.9011e+15 


3. 


.1623e-01 






3. 


.1004e+06 


3 


. 1623e- 


-07 




91 4.4p+1 R 


1 


.2284e+15 


5. 


.6234e-01 


4 


1 Q1 8p+flR 


1. 


.3546e+06 


5 


,6234e- 


-07 


5 


1 9nfip+1 4 


5. 


.2114e+14 


1, 


.0000e+00 




fifiRfip+DR 

. OOOUc T UJ 


5 


,8530e+05 


1 


.OOOOe- 


-06 


2 


. J. J. J.^ 


2 


.2168e+14 


1. 


,7783e+00 


g 


^34fip+D4 

. JOlUc T Ul 


2 


.4928e+05 


1 


,7783e- 


-06 


9 


1 n40p+1 3 


9 


.4594e+13 


3. 


.1623e+00 


2 




1 


.0423e+05 


3 


. 1623e- 


-06 


3 


838fip+1 3 


4 


,0507e+13 


5 


,6234e+00 


9 


41 73p+D3 


4. 


.2578e+04 


5 


.6234e- 


-06 






1 


,7409e+13 


1. 


.0000e+01 


3 


49Q1 p+H3 


1. 


.6897e+04 


1 


,0000e- 


-05 


6 


899Qp+1 9 


7 


.5091e+12 


1. 


.7783e+01 




9n4np+f)3 


6. 


.4693e+03 


1 


,7783e- 


-05 


2 


87R9p+1 9 


3. 


,2506e+12 


3. 


.1623e+01 


4 




2 


.3701e+03 


3 


. 1623e- 


-05 


1. 


2123e+12 


1. 


,4121e+12 


5. 


,6234e+01 


1 


.2859e+02 


8 


.2272e+02 


5 


.6234e- 


-05 


5 


1092e+ll 


6 


. 1560e+ll 


1. 


.0000e+02 


3. 


.8275e+01 


2 


.6741e+02 


1 


.OOOOe- 


-04 


2. 


.15280+11 


2 


,6929e+ll 


1 


,7783e+02 


1. 


.0511e+01 


8, 


,0233e+01 


1 


.7783e- 


-04 


9. 


0683e+10 


1. 


,1819e+ll 


3 


,1623e+02 


2. 


.6180e+00 


2. 


. 1844e+01 


3 


.1623e- 


-04 


3 


8186e+10 


5 


.2042e+10 


5. 


.6234e+02 


5. 


.79216-01 


5. 


,2861e+00 


5 


.6234e- 


-04 


1. 


6072e+10 


2 


.2985e+10 


1. 


.0000e+03 


1. 


.11016-01 


1. 


.1090e+00 


1 


.OOOOe- 


-03 


6. 


7607e+09 


1. 


.01816+10 


1, 


.7783e+03 


1 


,7883e-02 


1 


.9570e-01 


1 


.7783e- 


-03 


2 


8416e+09 


4 


.5213e+09 


3. 


.1623e+03 


2 


.3343e-03 


2. 


.8009e-02 


3 


. 1623e- 


-03 


1. 


, 1931e+09 


2 


.0124e+09 


5. 


.6234e+03 


2 


.3621e-04 


3 


.1107e-03 


5 


.6234e- 


-03 


5. 


0025e+08 


8. 


,9740e+08 


1. 


,0000e+04 


1. 


.7564e-05 


2 


.5414e-04 


1 


.OOOOe- 


-02 


2. 


0938e+08 


4 


,00706+08 


1. 


.7783e+04 


8. 


.9943e-07 


1. 


.4316e-05 


1 


,7783e- 


-02 


8. 


7436e+07 


1 


,7903e+08 


3. 


. 1623e+04 


2. 


.9327e-08 


5. 


.1409e-07 


3 


. 1623e- 


-02 


3. 


6404e+07 


7. 


,9973e+07 


5. 


.6234e+04 


5 


.5363e-10 


1 


.0702e-08 


5 


,6234e- 


-02 


1. 


5099e+07 


3. 


,5680e+07 


1. 


.0000e+05 


5. 


.3930e-12 


1. 


.15126-10 


1 


.OOOOe- 


-01 


6. 


. 2322e+06 


1 


.5878e+07 


1. 


.7783e+05 


2 


.3578e-14 


5 


,5658e-13 



TABLE B2 

Sample Gamma-Ray Spectra (a = 0) 



Flux 



Energy Inv. Comp. ee Brem. ep Brem. n° decay 
(eV) (photons s" 1 cm" 2 GeV- 1 ) 



1 


,0000e+06 


2 


.0062e- 


-10 


1 


,3845e- 


-10 


1 


.6562e- 


-10 


8 


.2454e- 


-18 


1. 


,7783e+06 


8 


.4541e- 


-11 


5 


.3800e- 


-11 


6 


.2697e- 


-11 


2 


.4513e- 


-17 


3 


. 1623e+06 


3 


,5616e- 


-11 


2 


.0632e- 


-11 


2 


.3393e- 


-11 


7 


.0967e- 


-17 


5 


.6234e+06 


1 


,5000e- 


-11 


7 


,7773e- 


-12 


8 


.5860e- 


-12 


1 


.9707e- 


-16 


1. 


,0000e+07 


6 


.3147e- 


-12 


2 


.8778e- 


-12 


3 


.10216- 


-12 


5 


. 1353e- 


-16 


1 


.7783e+07 


2 


.6568e- 


-12 


1 


,0460e- 


-12 


1 


. 1053e- 


-12 


1 


.2201e- 


-15 


3 


. 1623e+07 


1 


.11706- 


-12 


3 


,7410e- 


-13 


3 


,8916e- 


-13 


2 


,5548e- 


-15 


5 


,6234e+07 


4 


.6913e- 


-13 


1 


,3195e- 


-13 


1 


.3565e- 


-13 


2 


.7122e- 


-15 


1. 


,0000e+08 


1 


.9677e- 


-13 


4 


.6000e- 


-14 


4 


.6881e- 


-14 


2 


.7122e- 


-15 


1 


.7783e+08 


8 


.2389e- 


-14 


1 


.5880e- 


-14 


1 


.6085e- 


-14 


1 


.9824e- 


-15 


3 


, 1623e+08 


3 


.4416e- 


-14 


5 


.4375e- 


-15 


5 


.4842e- 


-15 


9 


.0019e- 


-16 


5 


,6234e+08 


1 


.4332e- 


-14 


1 


.8492e- 


-15 


1 


.8596e- 


-15 


3 


.6492e- 


-16 


1. 


.0000e+09 


5 


,9441e- 


-15 


6 


,2516e- 


-16 


6 


.2745e- 


-16 


1 


,3641e- 


-16 


1. 


.7783e+09 


2 


.4521e- 


-15 


2 


,1026e- 


-16 


2 


.1075e- 


-16 


4 


,8273e- 


-17 


3. 


, 1623e+09 


1 


.0045e- 


-15 


7 


,0378e- 


-17 


7 


.0481e- 


-17 


1. 


.6490e- 


-17 


5 


.6234e+09 


4 


.0771e- 


-16 


2 


.3445e- 


-17 


2 


.3468e- 


-17 


5 


.5070e- 


-18 


1 


.0000e+10 


1 


.6351e- 


-16 


7 


.7725e- 


-18 


7 


.7775e- 


-18 


1 


.81166- 


-18 


1. 


,7783e+10 


6 


.4567e- 


-17 


2 


.5626e- 


-18 


2 


.5637e- 


-18 


5 


.8906e- 


-19 


3 


. 1623e+10 


2 


.4988e- 


-17 


8 


.3928e- 


-19 


8 


.3953e- 


-19 


1. 


.8934e- 


-19 


5 


.6234e+10 


9 


.4223e- 


-18 


2 


,7251e- 


-19 


2 


.7257e- 


-19 


5 


.9973e- 


-20 


1. 


,0000e+ll 


3 


.4360e- 


-18 


8 


,7461e- 


-20 


8 


.7476e- 


-20 


1. 


.8582e- 


-20 


1. 


,7783e+ll 


1 


.2001e- 


-18 


2 


,7620e- 


-20 


2 


.7624e- 


-20 


5 


.5578e- 


-21 


3 


.16236+11 


3 


,9651e- 


-19 


8 


.5250e- 


-21 


8 


.5262e- 


-21 


1 


.5693e- 


-21 


5 


,6234e+ll 


1 


.21916- 


-19 


2 


. 5465e- 


-21 


2 


.5469e- 


-21 


4. 


.0316e- 


-22 


1. 


,0000e+12 


3 


.4130e- 


-20 


7 


.2548e- 


-22 


7 


.2560e- 


-22 


8 


.8745e- 


-23 


1 


.7783e+12 


8 


.4456e- 


-21 


1 


.9286e- 


■22 


1 


,9290e- 


-22 


1. 


.51616- 


-23 


3 


. 1623e+12 


1 


.7719e- 


-21 


4 


,6256e- 


-23 


4 


.6267e- 


-23 


1 


,7037e- 


-24 


5. 


,6234e+12 


2 


.9674e- 


-22 


9 


,4868e- 


-24 


9 


.4896e- 


-24 


9 


.4934e- 


-26 


1 


.0000e+13 


3 


.6207e- 


■23 


1 


.5226e- 


-24 


1 


.5232e- 


-24 


1 


.6037e- 


-27 


1 


.7783e+13 


2 


.7845e- 


-24 


1 


. 6433e- 


-25 


1 


.6443e- 


-25 


3 


.4422e- 


-30 


3 


,1623e+13 


1 


.0582e- 


-25 


9 


.1509e- 


-27 


9 


.1588e- 


-27 


1 


.9973e- 


-34 


5 


,6234e+13 


1 


.2918e- 


-27 


1 


. 6454e- 


-28 


1 


.6477e- 


-28 


1 


.9824e- 


-41 


1 


.0000e+14 


2 


,3039e- 


-30 


4 


. 1544e- 


-31 


4 


.1640e- 


-31 


2 


.4571e- 


-53 
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TABLE B3 
Sample Gamma-Ray Spectra (a = 0.05) 



Flux 

Energy Inv. Comp. ee Brem. ep Brem. n° decay 
(eV) (photons s" 1 cm" 2 GcV -1 ) 



1 


.0000e+06 


2 


,4116e- 


■10 


1. 


3856e- 


-10 


1 


,6572e- 


-10 


1 


. 1455e- 


-17 


1. 


,7783e+06 


1 


.0553e- 


-10 


5 


3856e- 


-11 


6 


.2750e- 


-11 


3 


. 1243e- 


-17 


3. 


, 1623e+06 


4 


.6340e- 


-11 


2. 


0662e- 


-11 


2 


.3422e- 


-11 


8 


.4245e- 


-17 


5 


.6234e+06 


2 


,0415e- 


-11 


7, 


7932e- 


-12 


8 


.6015e- 


-12 


2 


,2119e- 


-16 


1 


.0000e+07 


9 


,0225e- 


-12 


2. 


8864e- 


-12 


3 


. 1106e- 


-12 


5 


.5285e- 


-16 


1. 


,7783e+07 


3 


.9996e- 


-12 


1. 


0506e- 


-12 


1 


. 1099e- 


-12 


1 


.2763e- 


-15 


3 


. 1623e+07 


1 


.7779e- 


-12 


3. 


7663e- 


-13 


3 


.9168e- 


-13 


2 


.6248e- 


-15 


5 


.6234e+07 


7 


, 9236e- 


-13 


1. 


3332e- 


-13 


1 


.3701e- 


-13 


2 


.7832e- 


-15 


1. 


,0000e+08 


3 


. 5389e- 


-13 


4. 


. 6739e- 


-14 


4 


.7619e- 


-14 


2 


.7832e- 


-15 


1. 


,7783e+08 


1 


. 5832e- 


-13 


1 


6278e- 


-14 


1 


.6483e- 


-14 


2 


.0479e- 


-15 


3 


. 1623e+08 


7 


,0899e- 


-14 


5 


6505e- 


-15 


5 


.6971e- 


-15 


9 


.5028e- 


-16 


5 


.62346+08 


3 


. 1756e- 


-14 


1. 


9621e- 


-15 


1 


.9725e- 


-15 


3 


.9827e- 


-16 


1. 


.0000e+09 


1 


.4211e- 


-14 


6. 


,8424e- 


-16 


6 


.8653e- 


-16 


1. 


.5601e- 


-16 


1. 


.7783e+09 


6 


. 3444e- 


-15 


2. 


4051e- 


-16 


2 


.4101e- 


-16 


5 


.8712e- 


-17 


3 


. 1623e+09 


2 


,8206e- 


-15 


8. 


. 5278e- 


-17 


8 


,5385e- 


-17 


2 


. 1653e- 


-17 


5 


,6234e+09 


1 


. 2459e- 


-15 


3 


0488e- 


-17 


3 


.0513e- 


-17 


7 


.9260e- 


-18 


1 


.0000e+10 


5 


,4508e- 


-16 


1 


0983e- 


-17 


1 


.0988e- 


-17 


2 


.9006e- 


-18 


1 


.7783e+10 


2 


,3532e- 


-16 


3. 


.9805e- 


-18 


3 


.9818e- 


-18 


1 


.0644e- 


-18 


3 


.1623e+10 


9 


,9771e- 


-17 


1, 


4483e- 


-18 


1 


.4486e- 


-18 


3 


.9144e- 


-19 


5 


,6234e+10 


4 


. 1292e- 


-17 


5. 


.2743e- 


-19 


5 


.2752e- 


-19 


1. 


.4370e- 


-19 


1 


.0000e+ll 


1 


,6555e- 


-17 


1. 


.9143e- 


-19 


1 


.9145e- 


-19 


5 


.2220e- 


-20 


1 


.7783e+ll 


6 


.3677e- 


-18 


6. 


8852e- 


-20 


6 


.8859e- 


-20 


1 


.8516e- 


-20 


3 


, 1623e+ll 


2 


.3209e- 


-18 


2. 


4353e- 


-20 


2 


.4355e- 


-20 


6 


.2578e- 


-21 


5 


,6234e+ll 


7 


.8871e- 


-19 


8 


3813e- 


-21 


8 


.3822e- 


-21 


1 


.9414e- 


-21 


1 


.0000e+12 


2 


,4462e- 


-19 


2. 


7658e- 


-21 


2 


,7661e- 


-21 


5 


.2052e- 


-22 


1. 


,7783e+12 


6 


.7266e- 


-20 


8. 


5667e- 


-22 


8 


.5679e- 


-22 


1 


.0934e- 


-22 


3 


. 1623e+12 


1 


.5746e- 


-20 


2. 


4116e- 


-22 


2 


.4121e- 


-22 


1 


.5278e- 


-23 


5 


.6234e+12 


2 


,9588e- 


-21 


5. 


8625e- 


-23 


5 


.8639e- 


-23 


1 


.0730e- 


-24 


1. 


,0000e+13 


4 


.0826e- 


-22 


1 


. 1297e- 


-23 


1 


. 1301e- 


-23 


2 


.3213e- 


-26 


1 


,7783e+13 


3 


.5900e- 


-23 


1 


4881e- 


-24 


1 


.4888e- 


-24 


6 


.4948e- 


-29 


3 


. 1623e+13 


1 


. 5844e- 


-24 


1 


0307e- 


-25 


1 


.0315e- 


-25 


5 


.0034e- 


-33 


5 


,6234e+13 


2 


.2946e- 


-26 


2. 


3534e- 


-27 


2 


.3565e- 


-27 


6 


.7148e- 


-40 


1 


,0000e+14 


4 


.9900e- 


-29 


7. 


7046e- 


-30 


7 


.7217e- 


-30 


1 


. 1451e- 


-51 
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Note added in proof. — After this manuscript was submitted, papers by Kamae et al. (astro-ph/0605581) and by Kel- 
ner, Aharonian & Bugayov (astro-ph/0606058) appeared, presenting parameterized cross-sections for various particles 
produced of pp collisions. We plan to revise our pion-decay model to incorporate these new cross-sections. Details will 
be discussed in a subsequent paper. 



